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Abstract 

In classical tsunami-generation techniques, one neglects the dynamic sea bed dis- 
placement resulting from fracturing of a seismic fault. The present study takes into 
account these dynamic effects. Earth's crust is assumed to be a Kelvin- Voigt mate- 
rial. The seismic source is assumed to be a dislocation in a viscoelastic medium. The 
fluid motion is described by the classical nonlinear shallow water equations (NSWE) 
with time-dependent bathymetry. The viscoelastodynamic equations are solved by a 
finite-element method and the NSWE by a finite-volume scheme. A comparison be- 
tween static and dynamic tsunami-generation approaches is performed. The results 
of the numerical computations show differences between the two approaches and 
the dynamic effects could explain the complicated shapes of tsunami wave trains. 



1 Introduction 



The accuracy of the computation of the whole life of a tsunami, from generation to 
inundation, obviously depends on the construction of the initial condition. This is 
why the process of tsunami generation must be modelled as accurately as possible. 
Even though the constraint of being able to predict tsunami arrival time, height and 
location as fast as possible must be taken into account (in other words, a trade-off 
must be found between the precision and the speed of computation of the initial 
condition), we believe that so far the scientific community has not payed enough 
attention to the crucial subject of tsunami generation. 



After the pioneering work of Kajiura [10] it has become a common practice in the 
tsunami community to translate the static sea bed deformation generated by an un- 
derwater earthquake onto the free surface and let it propagate. We will refer to this 
method as passive approach. The validity of this technique was already discussed in 
[16,3]. Three-dimensional (3D) analytical expressions derived from Volterra's for- 
mula applied to the general study of dislocations [13,17] are used to construct the 
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static initial deformation. Similar analytical expressions for two-dimensional (2D) 
problems were also derived by Freund & Barnett [4] , who used the theory of analytic 
functions of a complex variable. The popularity of these analytic solutions can be 
explained by their relatively simple explicit form. Thus, their computation is easy 
and inexpensive. A feature of the solution of Freund & Barnett is that nonuniform 
slip distributions can be considered. In particular, slip distributions which remove 
the singular behavior of the internal stresses at the ends of the slip zone can be 
dealt with, simply by imposing the so-called smooth closure condition on the slip: 
the slip is zero at the ends. 

When simplifying hypotheses such as homogeneity or isotropy are removed, analyt- 
ical solutions are no longer available and the governing equations must be solved 
numerically. Static deformations caused by slip along a fault have been extensively 
simulated by Masterlark [14], who used several dislocation models based on the 
finite-element method (FEM) to estimate the importance of different physical hy- 
potheses. Anisotropy and heterogeneity turned out to be the most important factors 
in this type of modelling. Megna et al. [15] also used the FEM to compare numer- 
ical results with analytical solutions. However neither in [14] nor in [15] were the 
dynamical aspects and the coupling with hydrodynamics considered. Moreover the 
consequences for the resulting tsunami waves were not pointed out. 

When one uses as initial condition a static seismic source together with the trans- 
lation of the sea bed deformation onto the free surface, one neglects the rupture 
velocity and the rise time. Several studies have already been performed to under- 
stand wave formation due to different prescribed bottom motions by introducing 
either some type of rise time or some type of rupture velocity. For example, Todor- 
ovska & Trifunac [24] studied the generation of waves by a slowly spreading uplift 
of the bottom. Hammack [7] generated waves experimentally by raising or lowering 
a box at one end of a channel, and considered various laws for the rise or the fall of 
the box. In their review paper, Dutykh & Dias [2] generated waves theoretically by 
multiplying the static deformations caused by slip along a fault by various time laws: 
instantaneous, exponential, trigonometric and linear. Haskell [8] was one of the first 
to take into account the rupture velocity. In fact he considered both rise time, T, 
and rupture velocity, V. Consider the source shown in Figure 1. The two horizontal 
coordinates x and y, and the vertical coordinate z are denoted hy x — {x,y,z). Let 
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Fig. 1. Geometry of the source model. The fault has width W, length L, depth d and dip 
angle S. 
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b{x,t) denote the fault displacement function and bo{x) the final displacement. The 
following form for b{x, t) was considered by Haskell: 



b{x, t) 



t-c/v<o 

{bo/T){t-C/V)0<t-C/V <T 
bo t-C/V>T 



The coordinate C is a coordinate along the fault. Eq. (1) implies that ed, t — a, 
fracture front is established instantaneously over a width W of the y—axis at depth 
d. The front propagates unilaterally at constant velocity V over a length L cos 6 
of the a;— axis. At any fixed point on the fault plane the relative displacement 
increases at constant velocity from at i = C/^ ^ constant final value bo at 
t — T + C/V. Okumura & Kawata [19] used Haskell's approach to investigate the 
effects of rise time and rupture velocity on tsunami generation. They considered 
two cases of sea bottom motion: (i) with only rise time and (ii) with both rise 
time and rupture velocity. They found that the effects of rupture velocity are much 
smaller than those of rise time when the rise time is assumed to be long (over 
10 min). Ohmachi et al. [16] also considered rise time and rupture velocity but 
unfortunately the dynamics is not clearly explained in their paper. Apparently they 
did not solve the elastodynamic equations with the second-order time derivative (see 
next section). Another attempt to understand dynamical effects is that of Madariaga 
[12], who considers a dip-slip dislocation propagating in a half-space and solves the 
elastodynamic equations by using the double Laplace transform. The solution is 
elegant but rather complex. Unfortunately no plots of the deformation of the free 
surface are provided and the coupling with the water layer is not considered either. 
The present study can be considered as an attempt to understand the coupling 
between seismic faulting and hydrodynamics by integrating numerically the time- 
dependent elasticity equations as well as the time-dependent fiuid equations. The 
authors have already adrcsscd the problem of tsunami generation in [3,2]. The main 
feature of the present study is the use of a more realistic earthquake source model. 

The paper is organised as follows. In Section 2 we briefly describe the mathematical 
models, both for solid and fluid motions. Section 3 contains details on the numerical 
methods used to solve the governing equations. The numerical method for the solid 
motion is validated in Section 4. Section 5 provides a comparison between the 
traditional approach to tsunami generation (in which the static sea bed deformation 
is translated onto the free surface) and the more realistic approach of dynamic 
generation (in which the wavetrain is generated by the motion of the bottom). We 
reveal numerically that the dynamic generation of tsunamis can for example create 
a leading depression wave when one expects a leading elevation wave. 
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2 Mathematical models 



Even though the numerical results shown in this paper are for 2D configurations, 

the modeling is performed for 3D problems. The horizontal coordinates are denoted 
by X and y, while the vertical coordinate is denoted by z. The displacements are 
denoted by u^, Uy and u^- We use different origins along the vertical axis for the 
solid and fluid motions. In the earth domain, z — denotes the sea bed at rest 
(assumed to be flat) while in the fluid domain, z — denotes the sea surface at rest. 



2. 1 Dynamic fault model 

The fault is assumed to be inside a geological viscous medium. Earth's crust is 
assumed to be a viscoelastic material of density p. We choose the Kelvin- Voigt vis- 
cosity model [20] which consists in using complex elastic coefficients (with negative 
imaginary parts in order to dissipate wave energy). For isotropic media it means 
that the Lame coefficients have a nonpositive imaginary part: A* = — iAj, /x* = 
jij. — ijii, where A,.,//^ > and Ai,/ij > 0. The classical elasticity equations are 
obtained by choosing Aj = and /Xj = 0. On the time-scales relevant to our prob- 
lem, elasticity is sufficient and the assumption of a Kelvin- Voigt viscous material is 
unnecessary. But we keep it for the sake of completeness. 

Let cp and be the classical velocities for the propagation of P and S waves in a 
medium of density p: 



The factors Qp and Qs measure the viscosity of the geological medium. In this study 
we restrict our attention to the weakly viscous case (1/|(5p| -C 1 and l/IQsl ^ 1)- 

Let g, represent the stress tensor. The displacement field u{x,y,z,t) — {ux,Uy,Uz) 
satisfies the classical elastodynamic equations from continuum mechanics: 




Complex Lame coefficients yield complex velocities for wave propagation. 




where the coefficients Qp and Qs are defined as follows: 




V • £ = p 



(2) 
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It is common in seismology to assume that the stress tensor a is determined by 
Hooke's law through the strain tensor £ = | (Vm + V^wj . Therefore 

X*{V ■u)L + 2iJ,*g. (3) 

Thus, we come to the following linear viscoclastodynamic problem ^ : 

V • (A*(V ■ u)L + fJ.*{Vu + V'u)) = p-^. (4) 



Recall that the mechanical characteristics p, X* and p* can possibly depend on the 
spatial coordinates {x,y,z). However we will assume that they are constant in the 
numerical applications. 

The fault is modeled as a dislocation inside a viscoelastic material. This type of 
model is widely used for the interpretation of seismic motion. A dislocation is con- 
sidered as a surface (in 3D problems) or a line (in 2D problems) in a continuous 
medium where the displacement field is discontinuous. The displacement vector is 
increased by the amount of the Burgers vector b along any contour C enclosing the 
dislocation surface (or line), i.e. 

j) du = h. (5) 
c 



We let a dislocation run at speed V along a fault inclined at an angle 5 with respect 
to the horizontal. Rupture starts at position x — Q and z = —d (it is supposed to be 
infinitely long in the transverse direction) and propagates with constant rupture 
speed V for a finite time L/V in the direction 6 stopping at a distance L. Let ^ be a 
coordinate along the dislocation line. On the fault located in the interval < ( < L 
slip is assumed to be constant. The rise time is assumed to be 0. 



2.2 Fluid layer model 



Since the main purpose is to model tsunami generation processes and since tsunamis 
are long waves, it is natural to choose the nonlinear shallow water equations (NSWE) 
as hydrodynamic model. These equations are widely used in tsunami modelling, 
especially in codes for operational use [25,21]. The vahdity of the NSWE model and 
the question of the importance of dispersive effects have already been addressed by 
the authors in [11]. 

^ We use the prefix "visco-" due to the presence of the imaginary part in the Lame 
coefficients, which is responsible for small wave damping. 
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Let 7] denote the free-surface elevation with respect to the still water level z = 0, 
V = {vx,Vy) the horizontal velocity vector, g the acceleration due to gravity and 
z — —h{x, y, t) the bathymetry. The NSWE in dimensional form read 

I + ^ + f + 5V|.t + «v,^o. (6) 

The effect of the moving bottom appears in the source term —dh/dt in the first 
equation. The unknowns rj and v are functions of time and of the horizontal coordi- 
nates X and y. Since the NSWE are essentially obtained from depth-averaging the 
Euler equations, the dependence on the vertical coordinate z disappears from the 
equations. The coupling between the earth and fluid models is made through the 
function h{x, y, t) which describes the moving sea bottom bathymetry. 



3 Numerical methods 

In the present study we made two natural choices. The solid mechanics equations 
are solved using the FEM with fully implicit time integration, while for the hydrody- 
namic part we take advantage of the hyperbolic structure of the governing equations 
and use a solver based on the finite- volume scheme (see for example [1,11])- 

3. 1 Discretization of the viscoelastodynamic equations 

In order to apply the FEM one first rewrites the governing equation (4) in variational 
form. The time-derivative operator is discretized through a classical second-order 
finite-difference scheme. The method we use is fully implicit and has the advantage of 
being free of any Courant-Friedrichs-Lewy-type condition. In such problems implicit 
schemes become advantageous since the velocity of propagation of seismic waves is 
of the order of 3 to 4 km/s. We apply the P2 finite-element discretization procedure. 
For the numerical computations, the freely available code FreeFem-|— |- [9] is used. 

Let us say a few words about the boundary conditions and the treatment of the 
dislocation in the program. Concerning the boundary conditions, we assume that 
the sea bed is a free surface, that is ^ • n = at 2; = 0. The other boundaries 
are assumed to be fixed or, in other words, Dirichlet type boundary conditions 
M = are applied. The authors are aware of the reflective properties of this type 
of boundary conditions. In order to avoid the reflection of seismic waves along the 
boundaries during the simulation time, we take a computational domain which 
is sufficiently large. This approach is not computationally expensive since we use 
adaptive mesh algorithms [9] and in the regions far away from the fault, element 
sizes are considerably bigger than in the fault vicinity. 
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Next we discuss the implementation of the dislocation surface. Across the fault, the 
displacement field is discontinuous and satisfies the following relation: 



u^{x,t)—u {x,t) —b{x,t), (7) 



where the signs ± denote the upper and lower boundary of the dislocation surface, 
respectively. The propagation of Burger's vector along the fault is given by 

b{x,t)^boH{t-C/V), (8) 



where V is the rupture velocity, H the Heaviside unit step function and C a coordi- 
nate along the dislocation line. 



3.2 Finite-volume scheme for the nonlinear shallow-water equations 



Here we briefly describe the discretization of the model (6) by a standard cell- 
centered finite volume method. The system (6) can be written as 

^ + V-F(V) = 5, (9) 



where V = {r], u. v) is the vector of conservative variables (which coincide here with 
the physical variables), S = {—dh/dt, 0, 0) the source term, and, for every n G M^, 

F(V) -11= l^{h + ri){v-n), (^^\v\^ + gr]^ nj . (10) 



The Jacobian matrix A(V) • n is defined by 

A(V)-n-^^^^^^'''^-''^"'^^^'^^''^ 



dV 



(11) 



gn V <^n J 



The computational domain C is triangulated into a set of control volumes 
n = UkgtK. We integrate Eq. (9) on K: 

^[YdQ+ J2 [ F{y)-nKLda= [Sdfl, (12) 
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where nxi denotes the unit normal vector on K H L pointing into L and N{K) 
{LeT : area(i\: n L) 7^ 0} . Then, setting 



^ ^ K 



we approximate (12) by 

d^K , area(Lni^) 

where the numerical flux 

arealLriA J 

is explicitly computed by the FVCF formula of Ghidaglia et al. [5] : 

*(V, W;„) ^ F(V) . n + F(W) ■ n _ ^^^^^^^^^^ ^^^^ F(W) ■ n - F(V) ■ n ^^^^ 

Here the Jacobian matrix A„(yu) is defined in (11), ;u(V, W) is an arbitrary mean 
between V and W and sgn(M) is the matrix whose eigenvectors are those of M 
but whose eigenvalues are the signs of that of M. 

In our problem, the Jacobian matrix (11) has three distinct eigenvalues 

Ai3 = t/-n±Cs, A2 = 0, 



where Cg — g{h + 77) is the velocity of long gravity waves. The right and left 
eigenvectors can be easily computed. Then it is straightforward to compute the 
sign matrix appearing in (14) and the numerical scheme is thus completely defined. 

In this section we did not deal with boundary conditions. This is a complicated 
topic in finite- volume methods (see [6] for details) . 



4 Validation of the numerical method 



In this section we consider an analytic solution to the line dislocation problem in the 
static case. Use is made of the well-known result described for example by Freund 
& Barnett [4] or Okada [17]. In order to simplify the expressions, we only consider 
the 2D case (in other words, the fault is infinite in the y— direction). In fact the 
best expression is that given by equation (24) in [12]. We checked that it is in full 
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Fig. 2. Comparison between analytical and numerical solutions for a static 2D fault with 
a dip angle equal to 7r/2. 

agreement with the limit of Okada's solution as the width becomes infinite. The 
sketch of the domain is given in Figure 1. The fault has infinite width [W oo). 
Its length is L, its depth d and its dip angle S. 

In the present paper we only give the vertical displacement component along the 
free surface, since it plays the most important role in tsunami formation. It can be 
expressed as the difference between two contributions, that from a first dislocation 
located at the beginning of the fault and that from a second dislocation located at 
the end of the fault. Let dL = d — Lsm6. One has 



\bo\ 



X 



X 



L cos S 



dr 



(15) 



where 





1 








TT 



sin 5 arctan ^ 
d 



d{d cos 6 — xsin 6) 
x'^ + d^ 



(16) 



For the validation of our numerical method we chose a fault corresponding to a 
dip angle 6 = 7r/2. The values of the other parameters are given in Table 1. This 
problem was solved by FEM after neglecting the dynamic terms. The results of 
the comparison with solution (16) are given in Figure 2. Good qualitative and 
quantitative agreement can be seen. Megna et al. [15], who also considered static 
displacement due to uniform slip across a normal fault, compared the 2D FEM 
results with the analytical solution in the case of a normal fault. In their conclusion, 
they state that it is for the vertical component of the surface displacement that the 
discrepancies are the largest. 
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Earth free surface 




Fig. 3. Static deformation due to the dislocation corresponding to the parameters given 
in Table 1. 



Parameter 


Value 


Young modulus, E, GPa 


9.5 


Poisson ratio, ly 


0.27 


Damping coefficient, 


500 


Damping coefficient, fj,i 


200 


Fault depth, d, km 


4 


Fault length, L, km 


2 


Dip angle, 5, ° 


13 


Burger's vector length, 6o , m 


10 


Water depth (uniform), ho, m 


400 


Acceleration due to gravity, g, m/s^ 


9.8 



Table 1 

Parameters used in the numerical computations. The water depth and the spatial extent in 
the main direction of propagation were chosen so that dispersive effects can be neglected. 

5 Results of the simulation 



In this section, we use the set of physical parameters given in Table 1. The static 
sea bed deformation obtained with the analytical solution is depicted in Figure 3. 
Note that the only difference between Figures 2 and 3 is the value of the dip angle. 

In order to illustrate the numerical computations we chose several test cases of ac- 
tive/passive tsunami generation. The passive generation approach consists in trans- 
lating the static sea bed deformation onto the free surface and letting it propagate 
under gravity [10]. On the other hand, the active approach uses the bottom motion 
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for wave generation. We proceed by computing the first eight or fifteen seconds 
of the earthquake dynamics. Then the bottom configuration is assumed to remain 
frozen during the rest of the simulation. Concerning the dynamical aspects of rup- 
ture propagation, we consider the Heaviside-type approach (8) where the dislocation 
propagates along the fault with rupture velocity V. One could use instead a dis- 
location for which Burger's vector Bq is also space-dcpcndcnt. But the main goal 
of the present study is to make an attempt to include the dynamic displacement 
of the sea bed. In the dynamical approach, we consider three cases: the limiting 
case where the rupture velocity V is infinite, a fast event with V — 2.5 km/s and a 
slower event with V — 1 km/s. 

We show below the differences between the passive and the dynamic approaches. 
This question has already been addressed by the authors [3] in the framework of the 
linearized potential flow equations and of a simplifled model for bottom deformation. 

In the first comparison we use a strong coupling between the dynamic displacement 
of the sea bed and the fiuid layer equations and compare it with the passive ap- 
proach, in which the static solution shown in Figure 3 is translated onto the free 
surface as initial condition. The rupture velocity V is assumed to be infinite. More- 
over the earthquake dynamics is computed during the first eight seconds. The free 
surface at the beginning of the tsunami generation process is shown on Figure 4. 
Further steps of this process are given in Figures 5-6. One may have the impression 
that the passive solution does not evolve. In fact, the explanation lies in the pres- 
ence of two different time scales in this problem. The fast time scale is provided by 
the earthquake (P— and 5"— waves) and the slow one by water gravity waves. Since 
the active generation solution is directly coupled to the bottom dynamics, it evolves 
with the fast time scale. It is interesting to compare Figures 5 and 6. One can see 
that the active approach gives at the beginning an amplitude which is almost twice 
larger but the amplitudes become comparable a few seconds later. 

The free-surface elevations are computed until the wave enters the purely propa- 
gation stage. This corresponds to Figure 7. One notices that the resulting wave 
amplitude and velocity are almost the same. Of course the waveform is different. 
One can see as well that the location of the elevation wave is the same, while the 
depression wave is slightly shifted. It can be explained by the larger extent of the 
dynamic solution. Thus, we can conclude from this first comparison that if one is 
only interested in tsunami travel time or even in rough inundation zone estimation, 
the passive approach can be used. 

The second comparison focusses on the infiuence of the rupture velocity at two 
separate times (Figures 8 and 9). The differences between the fast and the relatively 
slow rupture velocities are small. 

The most interesting comparison is the third one, which focusses on the duration 
of the earthquake. Recall that our somewhat artificial definition of earthquake du- 
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Fig. 4. Water free surface at the beginning of the earthquake {t = 2s) according to two 
approaches of tsunami generation: passive versus active (with infinite rupture velocity). 













Fig. 5. Same as figure 4 for times t = As and t = 6s. 




Fig. 6. Same as figure 4 for times t = 7s and t = 10s. 
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Fig. 7. Same as figure 4 for times t = 150s and t = 250s. The wave is leaving the generation 
zone (left plot) and starting to propagate (right plot). 
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Fig. 8. Same as figure 7 (left plot) for two rupture velocities: V 
V = 2.5 km/s (right). 
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Fig. 9. Same as figure 8 at time t = 250s. 

ration is the time at which we stop the bottom motion. After that time, the sea 
bottom remains frozen. Figure 10 shows the effect of a longer earthquake. One sees 
that the shapes of the wave train obtained with the dynamic analysis look more 
complicated than those obtained with the passive analysis. In particular, the dis- 
tinction between leading elevation wave and leading depression wave is not as clear 
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Fig. 10. Water free surface at f. = 150s for a longer earthquake. The motion stops after 15s. 
The left plot compares the infinite rupture velocity solution with the passive approach. 
The right plot compares the case of a slow rupture {V = 1 km/s) with the passive case. 

when using the dynamical analysis. It could be an explanation for the discrepancies 
between modeled and recorded time series of water levels at various locations along 
the California coast for the 1960 Chilean tsunami. 



6 Conclusions 

An approach to model the dynamical character of sea bed deformations during an 
underwater earthquake was presented. The governing elastodynamic equations were 
solved by a finite-element method. The principal novelty of the present study is the 
coupling of the resulting displacement field with the hydrodynamic model. 

Two methods for tsunami generation have been compared: static versus dynamic. 
The computational results show that the dynamic approach leads to higher wa- 
ter levels in the near-fault area. These significant differences only occur during the 
first instants of the surface deformation and level off later on. However it was also 
observed that the shape of the wave train can be altered by dynamical effects. Con- 
sequently the distinction between leading elevation wave and leading depression 
wave may not be as clear as anticipated. Of course the present method is compu- 
tationally more expensive but there is an overall gain in accuracy. Not surprisingly 
more accurate tsunami computations require finer initial conditions such as those 
obtained by the active generation methodology used in the present study. 

In future work we intend to extend this modeling to three space dimensions since 
it is evident that the 2D computations presented in this article have little interest 
beyond academics. Moreover we intend to use the exact results of Tadepalli & 
Synolakis [22,23] to check how different the runup may be for the two different 
initial waves (resulting from the passive and dynamic seafloor displacements) in an 
idealized basin similar to the one used by Okal & Synolakis [18]. 
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